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We construct a mesoscale model of colloidal suspensions that contain solutes reversibly adsorbing 
onto the colloidal particle surfaces. The present model describes the coupled dynamics of the 
colloidal particles, the host fluid, and the solutes through the Newton-Euler equations of motion, the 
hydrodynamic equations, and the advection-diffusion equation, respectively. The solute adsorption 
is modeled through a square-well potential, which represents a short-range attractive interaction 
between a particle and a solute molecule. The present model is formulated to be solved through 
direct numerical simulations. Some numerical results are presented to validate the simulations. The 
present model enables investigations of solute adsorption effects in the presence of a fluid flow and 
an inhomogeneous solute concentration distribution. 


I. INTRODUCTION 

The dynamics of colloidal particles is appropriately in¬ 
vestigated by means of mesoscale modeling. Mesoscale 
models explicitly describe the motion of the colloidal 
particles and the host fluid, thereby enabling exact 
evaluation of the hydrodynamic interactions among the 
particles mediated by the motion of the surrounding 
fluid 0,0]. In many practical applications, solutes, such 
as electrolytes, surfactants, and adsorbing polymers, are 
added to control the dispersion stability of colloidal sus¬ 
pensions 01 - The mesoscale model can describe the 
influence of solutes on the dynamics of colloidal particles 
by additionally considering the temporal evolution of the 
solute concentration field. 

In the present study, we construct a mesoscale model 
of colloidal suspensions with an uncharged adsorbing so¬ 
lute such as a soluble polymer. The solute molecules 
reversibly adsorb onto the particle surfaces to form thick 
adsorption layers, whose overlapping generates interpar¬ 
ticle forces §5, 0. We describe the solute adsorption 
through a square-well potential, which represents a short- 
range attractive interaction between a particle and a so¬ 
lute molecule. This potential model reproduces the in¬ 
terparticle forces resulting from the overlapping of the 
adsorption layers, such as an attractive force induced by 
polymer bridging and a repulsive force induced by an 
increased osmotic pressure p-Q. A decrease in the so¬ 
lute concentration around the particles, namely negative 
adsorption, can be described by changing the particle- 
solute interaction to be repulsive, thereby generating an 
attractive interparticle force referred to as a depletion 
force [8]- The present model also reproduces a parti¬ 
cle motion caused by an inhomogeneity of solute con¬ 
centration distribution, such as diffusiophoresis 0,0- 
We describe the motion of the colloidal particles, the 
temporal evolution of the fluid flow field, and that of 
the solute concentration field through the Newton-Euler 
equations of motion, the hydrodynamic equations, the 


advection-diffusion equation, respectively. These gov¬ 
erning equations of the present model are coupled each 
other, an explicit form of which equations we derive. The 
solute transport is contributed by the fluid flow as ad- 
vective transport and is affected by the particles through 
boundary conditions on the diffusion flux derived from 
the solute impermeability and adsorptivity to the parti¬ 
cles. The coupled dynamics of the particles and the fluid 
is affected by the particle-solute adsorption interaction 
and the osmotic pressure in the adsorption layers. 

The present model is formulated to be solved through 
direct numerical simulations, of which scheme based on 
the immersed boundary method we have constructed 
[HB3. I n tbis method, the fluid flow field is also defined 
inside the particles, where a body force is added to satisfy 
no-slip boundary conditions at the particle surfaces. To 
indicate the domain occupied by the particles, we intro¬ 
duce a function that is equal to unity inside the particles 
and zero outside them. The two domains are smoothly 
connected through a thin interfacial region with finite 
thickness to represent the particle-fluid interface on the 
discrete computational grid points. In the present model, 
we also introduce another function to indicate the do¬ 
main occupied by the adsorption layers. The boundary 
conditions on the diffusion flux are imposed by use of the 
indicator functions for the particles and the adsorption 
layers. 

We examine the validity of the numerical simulations 
by comparing with some analytical solutions derived in 
this paper; one is a steady-state solute concentration dis¬ 
tribution around an isolated particle with imposing an 
external concentration gradient, and another is interpar¬ 
ticle force generated by the adsorption layer overlapping. 
We then investigate the influence of the solute adsorption 
on the steady-state sedimentation velocity of particles as 
an example of the situation where the dynamics of the 
particles, the fluid, and the solute molecules are coupled. 


II. MODEL 


We model a colloidal suspension as a system in which 
* tatsumi@sogo.t.u-tokyo.ac.jp rigid spherical particles with the same size and mass 
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FIG. 1. (Color online) Interaction potential between a particle 
and a solute molecule. The potential u is described by a thick 
solid line and consists of the hard-core exclusion potential 
u ex and the adsorption potential w ad , which potentials are 
described by broken lines. 


are dispersed in a Newtonian fluid. The fluid contains 
one solute species; the solute molecules reversibly adsorb 
onto the surface of the particles. The reversible adsorp¬ 
tion, namely physical adsorption, is described through 
a particle-solute interaction. For the sake of simplicity, 
we ignore solute-solute and solute-solvent interactions; 
in this assumption, the thermodynamic properties of the 
fluid is given by those of the ideal solution. 

In the present model, we consider the motion of the 
particles, the temporal evolution of the fluid flow field, 
and that of the solute concentration field. The motion 
of the i-th particle is described by the position Ri, the 
translational velocity V), and the angular velocity fl,. 
With the assumption that the fluid is incompressible, the 
fluid flow is described by the velocity field v(r,t) and the 
pressure field p(r, t), where r is the positional vector and 
t is the time. The solute concentration field c(r,t) is 
given as the number density of the solute molecules. 


A. Particle-solute interaction 


We describe the interaction between a particle and a 
solute molecule with a square-well potential, which is ex¬ 
pressed as follows: 


u(r) = u ex (r) + u ad (r), 

(1) 

f oo, 0 < r < a, 


(6X ( r ) = n 

(2) 

[ 0, r > a, 


{ -e, 0 < r < b, 


M = 

(3) 

1 0, r > b, 



where r represents the distance of the solute molecule 
from the center of the particle. A sketch of the square- 
well potential is shown in Fig. [L] The square-well po¬ 
tential u consists of two contributions: the hard-core ex¬ 
clusion potential u ex and the adsorption potential u ad . 
The hard-core exclusion potential u ex imposes the solute 
impermeability to the particle. The size of the solute 
molecules is assumed to be negligibly smaller than that of 
the particles, and hence the closest center-to-center dis¬ 
tance between a particle and a solute molecule is equal 
to the particle radius a. The square-well potential is 
characterized by the well-width w and the well-depth 
e. The adsorption layer is described as a spherical shell 
enclosed between two concentric spheres of radii a and 
b = a + w. The well-width and the well-depth correspond 
to the adsorption layer thickness and the adsorption en¬ 
ergy, respectively. Giving a negative value for the adsorp¬ 
tion energy e, negative adsorption can be described. We 
measure the adsorption energy by a non-dimensionalized 
value /3e, where /3 = (fc^T) -1 , Jcb is the Boltzmann con¬ 
stant, and T is the thermodynamic temperature. Assum¬ 
ing the additivity of potentials, the total potential acting 
on a solute molecule at the position r is given by 

U(r) = Y^ u ( r i), (4) 

i 

where r, = |rj|, and r,; = r —i?,. We also define the total 
adsorption potential as 

t/ ad (r) = 5> ad (r*)- (5) 


B. Indicator functions 


In the present formulation, fluid flow field is also de¬ 
fined inside the domain occupied by the particles. The 
domain occupied by the i-th particle Pi is given by 


Pi = {r | Xi(r,t;a) = 1}, 


where the function \i is a step function: 


Xi(r,t-,q) 


1, 0 < n < q, 
0, Ti > q. 


( 6 ) 

(7) 


We introduce a function <J>i to indicate the domain Pi as 


®i(r,t) = Xi{r,t;a). (8) 

A function indicating the union of all the domains {Pi} 
is also given by 

(9) 

i 

We assume that the particles interact through a hard¬ 
core exclusion potential to prevent them, namely the do¬ 
mains {Pi}, from overlapping one another, thereby lim¬ 
iting the range of the function $ to [0,1]. The particle 
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velocity field v p , which represents the rigid motion of the 
particles, is assigned to the velocity field inside the par¬ 
ticles: 


a virtual concentration field c*(r, t) as a continuous con¬ 
centration field from which the influence of the particles 
is eliminated: 


^ ^i x r i)- (10) 

i 

This formulation is equivalent to the imposition of no-slip 
boundary conditions on the velocity field at the particle 
surfaces. 

We also define the domain W) as one where the adsorp¬ 
tion potential is induced by the i -th particle, namely the 
domain inside the outer sphere composing the adsorption 
layer of the *-th particle: 

Wi = {r\xi(r,t-,b) = l}. (11) 


Since the domain Wi includes the domain Pi, the ad¬ 
sorption layer is represented by the domain Wj\Pi. We 
define a function 2* to indicate the domain Wi as the 
Boltzmann factor of the adsorption potential induced by 
the i -th particle: 


E»(r, t) = exp (-/3< d ), (12) 

= exp [pexi(r,t-,bj], (13) 


where uf d = u ad (ri). A function indicating the union of 
all the domains {Wi} is given by 


l(r,t) = exp(-/3t/ ad ), 


= exp 


P e 'Yhxi(r,t-,b) 




(14) 


(15) 


The domains {Wi} can overlap one another, and hence 
the function S is equal to exp(n/3e), where n is the over¬ 
lapping number of the domains {Wi}. The function 2 is 
unity outside the domains { W t } where n = 0. 

Under the ideal solution assumption, the chemical po¬ 
tential difference between the solute and solvent, n , is 
given as 

fi = kBTln— + U, (16) 

c* 

where c & is the solute concentration in the standard state. 
The concentration field in equilibrium with a fixed con¬ 
figuration of the particles is derived from the condition 
that /x is constant as 


c eq = c 0 exp(-/3t/), 

= (1 - ( t>)Sc 0 , 


(17) 

(18) 


where cq is the bulk solute concentration. The indicator 
functions give the influence of the particles on the concen¬ 
tration field. The solute concentration is zero inside the 
particles, namely {P;}, due to the solute impermeability 
to the particles and changes discontinuously by a factor of 
exp(n/3e) in the adsorption layers, namely {Wi\Pi}. For 
an arbitrary solute concentration field c(r,f), we define 


c = (1 — $)Sc*. (19) 

The virtual concentration field is also defined inside the 
particles. 


C. Advection-diffusion equation 

The solute concentration field obeys the following 
equation that describes the conservation of solute mass: 

^ + V • (cv s ) = 0, (20) 

where v s is the transfer velocity field of the solute. The 
diffusion flux J , which is defined as the relative mass 
flux to the barycentric fluid motion, is given in the linear 
order as [HI 

J = c{v s -v), (21) 

= -rV/i, (22) 

where T is the transport coefficient. Substituting Eq. © 

into Eq. (12^1) yields an expression of the diffusion flux as 

J = -D(l-$)EVc*, (23) 

where D = VksT/c is the diffusion coefficient of the so¬ 
lute; we assume a constant diffusion coefficient in the 
present model. Equation (E3l) indicates that the influ¬ 
ence of the particles on the diffusion flux appears in the 
same manner as that on the solute concentration field de¬ 
scribed by Eq. (fTUl) . This formula is equivalent to the im¬ 
position of boundary conditions, which are derived from 
the solute impermeability and adsorptivity to the parti¬ 
cles, on the diffusion flux. 

We consider the temporal evolution of the virtual con¬ 
centration field c* instead of the real concentration field 
c. Equation (l20l) is rewritten by substituting Eq. (fl9l) as 

^ [(1 - $)2c*] + V • [(1 - $)Sc*« s ] = 0. (24) 

Since there is no solute transport through the particle- 
fluid interface, mass conservation of the solute virtually 
distributed inside the particles is described as 

^(<J>Sc*) + V-(<l>5 C *i, p ) = 0. (25) 

Adding each side of Eqs. (THl) and E5l) together yields 
the following equation: 

|(Sc*) + V-(HcX) = 0, (26) 

where 


C V = C V 


■J*, 


(27) 
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J* = E~ 1 J = -D(l-$)Vc*. (28) 


D. Hydrodynamic equations 


Equation (1261) is rewritten as 

+ C * (W + ^' VH ) = °- (29) 

The indicator function Sj obeys the following advection 
equation: 

^ + V l -VE l = 0. (30) 

dt 

The advection velocity, which is not contributed by the 
angular velocity of the particle due to the spherical sym¬ 
metry of the adsorption potential u ad , is equal to the 
translational velocity of the particle. The second term 
in the left-hand side of Eq. (l29l) is rewritten by use of 
Eq. (13(1 as follows: 


The hydrodynamic equations of the incompressible 
fluid consist of the conservation equations of mass and 
momentum as 

V • v = 0, (34) 

P (^ +« ■ Vu) = V • (T + / s + /p, (35) 

where p is the mass density of the fluid. The stress tensor 
is given by 

<x = ~{p + 7r)J + r/[Vv + (V«) T ], (36) 

where rj is the shear viscosity. The osmotic pressure n, 
which is the pressure change in the adsorption layers ac¬ 
companied by the solute adsorption, is given by the van’t 
Hoff’s formula in the ideal solution assumption as 


dc* 

¥ + v '(^ 


— + V* • VS = V — ( 4- V* ■ V 

dt 8 ^ «• ' s 


<9S i 


dt 


= E§ (^-vo-vSi, 

i 

= -/3S^K^V))-V< d , 

i 

= -PeE^2(v* - Vi)-ri6(ri - b), (31) 


where fi = rj/rj. Substituting Eq. m into Eq. m 
yields the equation describing the temporal evolution of 
the virtual concentration field as 


dc* 

~dt 


V • (c*v*) 

jdec* ^2(v* - Vi) ■ fi5{n - b). 


(32) 


The right-hand side of Eq. (1321) represents that the in¬ 
flow/outflow of the virtual concentration to/from the do¬ 
main , namely the adsorption/desorption of the so¬ 
lute, results in the dissipation/production of the vir¬ 
tual concentration at the surface dWi. This dissipa¬ 
tion/production occurs to satisfy mass conservation of 
the solute, whose real concentration in the adsorption 
layer is given by the multiplication of the virtual concen¬ 
tration with a factor of exp(n/3e). The transport veloc¬ 
ity of the virtual concentration through the surface dW t , 
which moves with the velocity V,, is given by v* — V. 

We consider thermal fluctuations, which affect the 
mesoscale dynamics of colloidal particles, through a ran¬ 
dom diffusion flux £ added to Eq. m- The random 
diffusion flux is a vector-valued stochastic variable that 
satisfies the fluctuation-dissipation relation fl4j |: 


{Ca{r,t)Cp(r',t')) = 2 Dc5 a p6(r - r')6{t - t'). (33) 


Tr = k B T(E-l)c*. (37) 

The body force fp is added to constrain the velocity field 
inside the particles to be equal to v p given in Eq. m- 
The body force fs is exerted by the particle-solute ad¬ 
sorption interaction and is expressed as 

fs = — cVC/ ad , (38) 

= k B T{ l-$)c*VS. (39) 

In Eq. 63, the effects of the solute adsorption appear 
through the chemical potential given by Eq. (fl6l) . of 
which the first term (the entropic contribution) and the 
second term (the energetic contribution) yield the os¬ 
motic pressure gradient — V7r and the body force fs, 
respectively. 

Thermal fluctuations are also considered for the ve¬ 
locity field through a random stress tensor s added to 
Eq. m- The random stress tensor is a tensor-valued 
stochastic variable satisfying the fluctuation-dissipation 
relation fl4| : 

( s a p(r,t)s IJ ,„(r',t')) 

= 2k B Tr)(8 ailj 8p v + 8 a ^5 0IJ ,)S(r - r')S(t - t’). (40) 


E. Newton-Euler equations of motion 

The motion of the colloidal particles is governed by the 
Newton-Euler equations of motion: 

M ~V i = F l H + F l c + F f, ^Ri = Vi, (41) 

h • = IVf. (42) 

The i-th particle has the mass M; and the inertia tensor 
Jj. The force Ff is exerted by direct interactions among 
the particles and prevents the overlap of the particles. 
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The hydrodynamic force F t H and torque N/P stem from 
the momentum exchange between the particle and fluid: 

F? = [ cr-dS, (43) 

J a Pi 


N? = nx(<T- d S), (44) 

■Id Pi 

where d5 is the area element vector. The integrations are 
performed over the surface of the particle dPi. The force 
Ff is the net force through the adsorption interaction 
between the particle and solute molecules: 

Ff = j cVufdr, (45) 

= k B T(e f3s - 1) [ #-(l-$)c*d S, (46) 

JdWi 

which describes the integration of the stress resulting 
from the discontinuous solute concentration change at 
the surface of the adsorption layer d W). Since the direc¬ 
tion of each force exerted by solute molecules, Vuf d , is 
normal to the particle surface, no torque is induced on 
the particle by the particle-solute adsorption interaction. 


III. ANALYTICAL SOLUTIONS 
A. Solute distribution disturbed by a particle 

Suppose that a constant solute concentration gradient 
K is imposed on a quiescent fluid, thereby generating 
a linear concentration distribution given by c = cq + 
K r in the absence of particles. We derive a disturbed 
concentration distribution in steady state when a single 
particle is fixed in the fluid. The origin of the coordinate 
system is set at the center of the particle, and hence the 
unit vector r = r/r with r = |r| is normal to the surfaces 
of the particle and the adsorption layer. In the present 
situation, a steady-state concentration field obeys 


V • J = 0, 

(47) 

with the boundary condition 


c* — > cq + K ■ r as r — > 00 . 

(48) 

Equation (l47l) is reformulated to the Laplace’s equation 
with the boundary conditions at the surfaces of the par¬ 
ticle and the adsorption layer: 

vV = 0, 

(49) 

(f • Vc*) r=a+ = 0, 

(50) 

(f • Vc*) r= 6+ = e /?£ (f • Vc*) r=b ~, 

(51) 

(r x Vc*) r=b + =(rx Vc*) r=6 -. 

(52) 


The conditions given by Eqs. 0 and m are derived 
from the discontinuity of the diffusion flux described in 
Eq. (E^l) . The condition given by Eq. (1521) . which is de¬ 
rived from the irrotatinality of Vc*, represents the con¬ 
tinuity of the concentration gradient in the tangential 
directions to the surface of the adsorption layer. The 
solution of Eqs. EHHED is obtained as 


c - c 0 = 


a ' 1 + 2^3 ) ( K ' r )’ a < r < bl 


b 3 \ 

1 + 'y^)( K -r), r>b , 


2 r 3 J 


(53) 


where 


3 _ 2 (-A + B) 

01 ~ A + 2B’ 7 ““ A + 2B 1 


(54) 


W '( 1 -S)' B = 1 + w- (55) 

The solutions with the adsorption layer thickness 
w/a = 0.5 and some values of the adsorption energy /3e 
are shown in Fig. [4] The result of non-adsorbing solute 
(/?£ = 0) indicates that the solute diffusion from high to 
low concentration is hampered by the particle, thereby 
increasing and decreasing the solute concentration in 
front of and behind the particle in the direction of K, 
respectively. The disturbance is diminished/enhanced by 
the positive/negative solute adsorption. This change re¬ 
flects the influence of the solute adsorption on the dif¬ 
fusion flux around the particle described by Eq. (l23l) . 
which indicates that the positive/negative adsorption in¬ 
creases/decreases the magnitude of the diffusion flux in 
the adsorption layer by a factor of exp(/3e). 


B. Interparticle force generated by adsorption 
layer ovelapping 

Suppose two particles arranged in a fluid as described 
in Fig. [2] The center-to-center vector from the particles 
1 to 2 is given by R = R 2 — Ri , the unit vector of which 
direction is expressed as R = R/R with R = |fi|. We 
consider an equilibrium solute concentration distribution 
with fixing the positions of the particles. When each 
particle is isolated from the other, zero net force acts on 
each particle because of the spherically symmetric solute 
concentration distribution around the particle. On the 
other hand, when the adsorption layers of the particles 
overlap each other, the spherical symmetry of the solute 
concentration distribution is broken, and hence non-zero 
net force acts on each particle. The forces acting on 
the particle 2, namely the force induced by the particle- 
solute adsorption interaction F.f and that induced by 
the osmotic pressure F.j^, are respectively evaluated by 
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(a) 



Particle 1 Particle 2 


FIG. 2. Configuration of two particles. The particles are sit¬ 
uated in a fluid with a center-to-center distance R = |i?| ■ (a) 
The adsorption layers overlap each other without one particle 
facing to the adsorption layer of the other as a + b < R < 2b. 
(b) One particle faces to the adsorption layer of the other as 
R < a + b. 



FIG. 3. (Color online) Force acting on a particle, whose ad¬ 
sorption layer overlaps with that of another particle, as a func¬ 
tion of the center-to-center distance between the particles for 
various adsorption energies fie . The adsorption layer thick¬ 
ness is set as w/a = 0.5. The negative and positive values 
correspond to attractive and repulsive forces, respectively. 


Eqs. (H7>1) and (T£J1) as 


The net force acting on the particle 2 is given by 


F 2 s = k B Tc 0 (e? £ - 1) / Sr(l - $)dS, 

JdW 2 

= -k B Tc o (e 0£ - 1) 

x [e 0e (Sww ~ Spw ) — Sww] R , (56) 

F 2 h = - f ndS , (57) 

JdP 2 

= -k B Tc o e 0£ [ Hid S, 

JdP 2 

= k B Tc o {e 20£ - e pe )S PW R , (58) 

where Sww and Spw are the circular areas enclosed by 
the edges of the surfaces dW 2 T W\ and dP 2 IT W\ (or 
dW 2 T Pi), respectively. In the derivation of Eqs. m 
and (1581) , we employ the following relations derived from 
the axisymmetry about the R axis: 


/ d5 = SwwRi 
Jaw 2 r\w 1 


(59) 


f dS= [ d S = S PW R- (60) 

JdP 2 rw i Jdw 2 np i 

The areas Sww and Spw are expressed as functions of 
the center-to-center distance between the particles R : 

Sww = — R 2 ) as R < 2b, (61) 

Spw = ^[^b 2 R 2 -(R 2 + b 2 ~a 2 ) 2 ] 

as R < a + b. (62) 


F 2 = Ff + F 2 S 
= k B Tc 0 (e^ s - 1) 

x [2 e^ £ Spw — (e^ — 1 )<S'vixvk] R-, (63) 

which is the same as the result obtained in Ref. [f|. 

Figure [5] shows the force acting on the particle 2 for 
the adsorption energies /3e = ±0.25, ±0.5 with the ad¬ 
sorption layer thickness w/a = 0.5. We first consider the 
positive adsorption described by fie > 0. The overlap¬ 
ping of the adsorption layers as R < 2b (where 2b = 3) 
leads to an increased solute concentration at the over¬ 
lapped adsorption layer, toward which the particles are 
attracted through the particle-solute adsorption interac¬ 
tion (Fig. (2Jr). In other words, an attractive force given 
by F 2 acts on the particle, which force describes an at¬ 
tractive force induced by polymer bridging. When two 
particles get closer as R < a ± b (where a ± b = 2.5), one 
particle faces to the adsorption layer of the other and 
goes through a high osmotic pressure at the overlapped 
adsorption layer (Fig. [2jo) , thereby imposing a repulsive 
force given by F 2 on the particle. Consequently, there 
exists an equilibrium interparticle distance as R/a ~ 2.43 
and 2.35 for /3e = 0.25 and 0.5, respectively; the particles 
will assemble with some interparticle separation. On the 
other hand, when we consider the negative adsorption de¬ 
scribed by fie < 0, a decreased solute concentration at the 
overlapped negative adsorption layers, where both of the 
particle-solute repulsive force and the osmotic pressure 
is diminished, imposes a totally attractive interparticle 
force. This attractive force is referred to as a depletion 
force 
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IV. NUMERICAL RESULTS 


Direct numerical simulations of the present model are 
performed by use of the scheme given in our previous 
papers [D], Q5 • To represent the indicator functions on 
the discrete computational grids, we use a continuous 
function as a substitute for the step function given by 
Eq. © [3: 


Xifaq) = 


/[(g + £/2) - n\ 


f[(q + £/2) - n] + f[n - (q - £/2)] 


,(64) 


f(x) 


exp(—A 2 /;r 2 ), x > 0, 
0, x < 0, 


(65) 


where A is the spacing of simulation grids. The func¬ 
tion given by Eq. (l64l) represents that the domains inside 
(d> = 1) and outside (<f> = 0) the particles are smoothly 
connected by a thin interface of thickness £, which is set 
as £/A = 2 in all of the simulations in the present study. 
In this description, the particle-fluid interface has a finite 
volume and thereby is supported by multiple grid points. 

The equations of motion of particles given by Eqs. m 
and (1121) are solved with the Euler scheme; time evolution 
of the particle position is solved with the Crank-Nicolson 
scheme. The hydrodynamic equations given by Eqs. AMD 
and (1551) are solved using the SIMPLEST (semi-implicit 
method for pressure-linked equations shortened) [L6|. 
The detailed expression of the constraint body force f p 
is described in the previous paper [T3 | . To solve the 
diffusion-advection equation of the virtual concentration 
field given by Eq. (1521) . the flux term in the left-hand side 
is rewritten by use of the solenoidal condition Eq. (1521) 
as 


V • (c*v*) = v ■ Vc* + V • J*. (66) 

The first and second terms in the right-hand side repre¬ 
sent the advection and the diffusion of the solute, respec¬ 
tively. The discretization in space of Eq. (1521) is carried 
out; a second-order central differencing is employd for the 
diffusion term, while a third-order upwind differencing is 
employed for the advection term and the right-hand side 
of the equation. The time integration is performed by the 
second-order Runge-Kutta scheme (the Heun scheme). 
In the following simulations, we neglect the thermal fluc¬ 
tuations introduced through Eqs. (1551) and (1201) . 


A. Solute distribution disturbed by a particle 

To validate the numerical simulation of the advection- 
diffusion equation given by Eq. (1521) . we calculated a so¬ 
lute concentration distribution around a particle with im¬ 
posing a concentration gradient. The analytical solution 
is derived in Hill Al The simulation box is set as a cube 
of side length L = 10a. A single particle is fixed on the 
center of the simulation box. We introduce a Cartesian 



x la 
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FIG. 4. (Color online) Virtual solute concentration field 
c* around a particle with the imposed concentration gradi¬ 
ent K = — 0.1(co/a)i. The adsorption layer thickness is 
w/a = 0.5. (a) Contour lines of the concentration field at 
the plane of z = 0 for fie = 0.5. The color lines represent the 
analytical solution, while the black dashed-two dotted lines 
represent the simulation result. The surfaces of the particle 
and the adsorption layer are described by thick and thin bro¬ 
ken lines, respectively, (b) Concentration distribution on the 
line of y = z = 0. The adsorption energy takes the values 
of fie = 0, ±0.5. The solid lines represent analytical solu¬ 
tions. The dashed double-dotted line represents the concen¬ 
tration distribution in the absence of the particle, c* = |DC|a;. 
The thick and thin broken lines represent the positions where 
surfaces of the particle and the adsorption layer are present, 
respectively. 


coordinate system with the origin at the center of the 
simulation box and the axes parallel to the sides of the 
simulation box. The imposed concentration gradient is 
K = —0.1(co/a)i, where x is the unit vector of x di¬ 
rection. The concentration at the boundary surfaces of 
x = ±L/2 is fixed as Co =F \K\L/2; periodic boundary 
conditions are imposed on the y and z directions. 
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FIG. 5. (Color online) Simulation estimation of force acting 
on a particle induced by the adsorption layer overlapping with 
various grid resolutions of the particle a/ A. The adsorption 
layer thickness and adsorption energy are w/a = 0.5 and fie = 
0.5, respectively. 


The virtual concentration field c* around the particle 
in steady state is shown in Fig. [I] The adsorption layer 
thickness and the adsorption energy take the values of 
w/a = 0.5 and fie = 0,±0.5, respectively. Excellent 
agreement with the analytical solutions is obtained. In¬ 
fluence from the periodic boundary conditions is slight 
because the order ratio of the disturbance caused by the 
particle to the non-disturbed concentration field is pro¬ 
portional to r -3 as described by Eq. m- 



FIG. 6. (Color online) Sedimentation velocity of a particle 
with different adsorption energies fie as a function of the ad¬ 
sorption layer thickness w/a. The velocity is normalized by 
the sedimentation velocity with no solute adsorption (fie = 0), 
which is given by Vo. The adsorption energy takes the values 
of /5e = ±0.25, ±0.5. 
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B. Interparticle force generated by adsorption 
layer ovelapping 

To validate the numerical evaluation of force induced 
on the particles by the solute adsorption, we calculated 
the interparticle force between two particles generated by 
the adsorption layer overlapping. The analytical solution 
is derived in 11IIIBI Figure [5] shows the simulation results 
with various grid resolutions of the particle a/ A. The 
adsorption layer thickness and the adsorption energy are 
set as w/a = 0.5 and fie = 0.5, respectively. A remark¬ 
able deviation from the analytical solution is observed 
around the kink at the minimum point. This discrepancy 
is attributed to the smoothed particle-fluid interface in 
the present simulation. An increase in the grid resolu¬ 
tion a /A leads to the reduction in the deviation of the 
simulation result from the analytical solution. 


C. Solute adsorption effects in sedimentation 

We present the simulation results of a steady-state sed¬ 
imentation velocity of particles arranged in a periodic ar- 
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FIG. 7. (Color online) Deviation of virtual solute concentra¬ 
tion c* — Co and velocity field v around a particle, with the 
adsorption layer thickness w/a = 0.5. The deviation of the 
virtual concentration is scaled by the bulk concentration co- 
The adsorption energy takes the values of (a) fie = 0.5 and (b) 
fie = —0.5. The cross-sections presented here are parallel to 
the particle motion direction and include the particle center. 
The particle is represented by a black circle. The direction of 
the particle velocity is to the downward in the pictures. The 
deviation of the virtual solute concentration is described by 
a color scale, which goes from negative (darker) to positive 
(lighter) deviation. 
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ray with an adsorbing solute; in this situation, coupled 
dynamics of the particles, the fluid, and the solute is to 
be considered. The simulation box is set as a cube of 
side length L = 8a with periodic boundary conditions, 
thereby forming a periodic array of the particles. In 
the reduced units as A = p = rj = 1, the bulk solute 
concentration is set as cofcsT = 0.0244, and the diffu¬ 
sion coefficient is set as D = 0.1 (the Schmidt number 
is Sc = p/pD = 10). The particles are continuously 
subjected to a constant force, of which magnitude is set 
such that the particle Reynolds number with no solute 
adsorption [fie = 0) is Re p = paVo/r] = 8 x 10~ 1 2 3 4 5 6 7 8 9 10 . 

The steady-state sedimentation velocity of the parti¬ 
cles as a function of the adsorption layer thickness is 
shown in Fig. [G] The adsorption energy takes the val¬ 
ues of fie = ±0.25, ±0.5. The sedimentation velocity is 
reduced with the increase in the adsorption layer thick¬ 
ness for both positive and negative adsorption, and hence 
the particle Reynolds number is less than 8 x 10 -3 in 
all the cases. The virtual concentration distribution and 
velocity field around one of the settling particles with 
w/a = 0.5 are depicted in Fig. [7] The positive adsorp¬ 
tion described by fie = 0.5 (Fig. (TJr) leads to a decreased 
and increased virtual concentration in front of and be¬ 
hind the particle in the direction of the particle motion, 
respectively, while the negative adsorption described by 
fie = —0.5 (Fig. [7Jr) leads to an opposite change in the 
virtual concentration. In both cases, the generated asym¬ 
metric concentration distribution induces a force on the 
particle in the direction opposite to the particle motion, 
thereby reducing the sedimentation velocity. The larger 
force induced on the particle according to the thicker ad¬ 
sorption layer is due to the larger surface area of the ad¬ 
sorption layer where solute molecules exert force on the 
particle. The generation of the asymmetric concentra¬ 
tion distribution can be explained in terms of Eq. (1321) : 


the solute transported slowly relative to the particle mo¬ 
tion is adsorbed and desorbed in front of and behind the 
particle. Similar phenomenon has been investigated for 
charged particles in an electrolyte solution 00- 

V. CONCLUSION 

We constructed a mesoscale model of colloidal suspen¬ 
sions containing an adsorbing solute. We derived equa¬ 
tions describing the coupled dynamics of the colloidal 
particles, the host fluid, and the solute molecules; these 
equations were formulated to be solved through direct 
numerical simulations. We performed the validation of 
the simulations by comparing with the analytical solu¬ 
tions derived in the present paper for two situations: a 
steady-state solute concentration distribution around a 
particle with imposing an external concentration gradi¬ 
ent and interparticle force generated by the adsorption 
layer overlapping. We then computed the steady-state 
sedimentation velocity of particles arranged in a periodic 
array. The reduction of the sedimentation velocity was 
observed according to an asymmetric solute concentra¬ 
tion distribution around the particles, which asymmetric 
distribution was generated by the solute adsorption and 
desorption at the surfaces of the adsorption layers. 

The present model is appropriate to investigate solute 
adsorption effects on the dynamics of colloidal particles 
in the presence of a fluid flow and an inhomogeneous so¬ 
lute concentration distribution. As an interesting appli¬ 
cation of the present model, we will be able to investigate 
the structure formation of colloidal particles in a drying 
process of colloidal suspensions containing adsorbing so¬ 
lutes. In this case, the present model is required to be 
developed to consider the gas-liquid free surface motion 
according to evaporation [19J, beneath which surface a 
high concentration layer of solutes can appear. 
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